<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Design a 1/f spectrum shaping (pink-noise) filter</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/filter_design/html/one_over_f_filter.html">
<link rel="stylesheet" href="../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Design a 1/f spectrum shaping (pink-noise) filter</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% "Filter design" lecture notes (EE364) by S. Boyd</span>
<span class="comment">% "FIR filter design via spectral factorization and convex optimization"</span>
<span class="comment">% by S.-P. Wu, S. Boyd, and L. Vandenberghe</span>
<span class="comment">% (a figure is generated)</span>
<span class="comment">%</span>
<span class="comment">% Designs a log-Chebychev filter magnitude design given as:</span>
<span class="comment">%</span>
<span class="comment">%   minimize   max| log|H(w)| - log D(w) |   for w in [0,pi]</span>
<span class="comment">%</span>
<span class="comment">% where variables are impulse response coefficients h, and data</span>
<span class="comment">% is the desired frequency response magnitude D(w).</span>
<span class="comment">%</span>
<span class="comment">% We can express and solve the log-Chebychev problem above as</span>
<span class="comment">%</span>
<span class="comment">%   minimize   max( R(w)/D(w)^2, D(w)^2/R(w) )</span>
<span class="comment">%       s.t.   R(w) = |H(w)|^2   for w in [0,pi]</span>
<span class="comment">%</span>
<span class="comment">% where we now use the auto-correlation coeffients r as variables.</span>
<span class="comment">%</span>
<span class="comment">% As an example we consider the 1/sqrt(w) spectrum shaping filter</span>
<span class="comment">% (the so-called pink-noise filter) where D(w) = 1/sqrt(w).</span>
<span class="comment">% Here we use a logarithmically sampled freq range w = [0.01*pi,pi].</span>
<span class="comment">%</span>
<span class="comment">% Written for CVX by Almir Mutapcic 02/02/06</span>

<span class="comment">% parameters</span>
n = 40;      <span class="comment">% filter order</span>
m = 15*n;    <span class="comment">% frequency discretization (rule-of-thumb)</span>

<span class="comment">% log-space frequency specification</span>
wa = 0.01*pi; wb = pi;
wl = logspace(log10(wa),log10(wb),m)';

<span class="comment">% desired frequency response (pink-noise filter)</span>
D = 1./sqrt(wl);

<span class="comment">% matrix of cosines to compute the power spectrum</span>
Al = [ones(m,1) 2*cos(kron(wl,[1:n-1]))];

<span class="comment">% solve the problem using cvx</span>
cvx_begin
  variable <span class="string">r(n,1)</span>   <span class="comment">% auto-correlation coefficients</span>
  variable <span class="string">R(m,1)</span>   <span class="comment">% power spectrum</span>

  <span class="comment">% log-chebychev minimax design</span>
  minimize( max( max( [R./(D.^2)  (D.^2).*inv_pos(R)]' ) ) )
  subject <span class="string">to</span>
     <span class="comment">% power spectrum constraint</span>
     R == Al*r;
cvx_end

<span class="comment">% check if problem was successfully solved</span>
disp([<span class="string">'Problem is '</span> cvx_status])
<span class="keyword">if</span> ~strfind(cvx_status,<span class="string">'Solved'</span>)
  <span class="keyword">return</span>
<span class="keyword">end</span>

<span class="comment">% spectral factorization</span>
h = spectral_fact(r);

<span class="comment">% figures</span>
figure(1)
H = exp(-j*kron(wl,[0:n-1]))*h;
loglog(wl,abs(H),wl,D,<span class="string">'r--'</span>)
set(gca,<span class="string">'XLim'</span>,[wa pi])
xlabel(<span class="string">'freq w'</span>)
ylabel(<span class="string">'mag H(w) and D(w)'</span>)
legend(<span class="string">'optimized'</span>,<span class="string">'desired'</span>)
</pre>
<a id="output"></a>
<pre class="codeoutput">
 
Calling sedumi: 3641 variables, 2400 equality constraints
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
Split 41 free variables
eqs m = 2400, order n = 3083, dim = 4283, blocks = 601
nnz(A) = 5400 + 49200, nnz(ADA) = 9600, nnz(L) = 6000
Handling 82 + 0 dense columns.
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.04E-01 0.000
  1 :   3.49E+01 8.25E-02 0.000 0.2042 0.9000 0.9000  -0.98  1  1  6.4E+01
  2 :   2.06E+02 7.52E-03 0.000 0.0911 0.9900 0.9900  -0.82  1  1  2.8E+01
  3 :   9.16E+01 3.21E-03 0.000 0.4272 0.9000 0.9000   1.69  1  1  7.7E+00
  4 :   3.04E+01 9.91E-04 0.000 0.3088 0.9000 0.9000   3.48  1  1  8.0E-01
  5 :   8.04E+00 4.32E-04 0.000 0.4361 0.9000 0.9000   4.95  1  1  1.0E-01
  6 :   3.44E+00 2.23E-04 0.000 0.5151 0.9000 0.9000   2.91  1  1  3.2E-02
  7 :   2.17E+00 1.01E-04 0.000 0.4548 0.9000 0.9000   1.53  1  1  1.3E-02
  8 :   1.68E+00 5.00E-05 0.000 0.4938 0.9000 0.9000   1.08  1  1  6.4E-03
  9 :   1.45E+00 3.11E-05 0.000 0.6218 0.9000 0.9000   0.92  1  1  4.1E-03
 10 :   1.27E+00 1.52E-05 0.000 0.4900 0.9000 0.9000   0.91  1  1  2.1E-03
 11 :   1.23E+00 8.33E-06 0.000 0.5466 0.9000 0.9000   1.02  1  1  1.1E-03
 12 :   1.20E+00 2.57E-06 0.000 0.3089 0.9016 0.9000   1.00  1  1  3.9E-04
 13 :   1.19E+00 7.43E-07 0.000 0.2887 0.9044 0.9000   1.00  1  1  1.2E-04
 14 :   1.19E+00 2.24E-07 0.000 0.3009 0.9012 0.9000   1.00  1  1  3.8E-05
 15 :   1.19E+00 2.61E-08 0.000 0.1168 0.9000 0.0000   1.00  1  1  1.3E-05
 16 :   1.19E+00 4.21E-09 0.000 0.1611 0.9292 0.9000   1.00  5  5  3.5E-06
 17 :   1.19E+00 2.24E-10 0.176 0.0532 0.9900 0.8463   1.00  6  6  3.3E-07
 18 :   1.19E+00 3.95E-12 0.000 0.0176 0.9900 0.9900   1.00 14 14  6.2E-09

iter seconds digits       c*x               b*y
 18      8.5   Inf  1.1873313477e+00  1.1873324283e+00
|Ax-b| =   2.5e-08, [Ay-c]_+ =   2.1E-09, |x|=  2.7e+02, |y|=  1.1e+00

Detailed timing (sec)
   Pre          IPM          Post
4.000E-02    8.540E+00    3.000E-02    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=7, |skip| = 2, ||L.L|| = 3.72784.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.18733
Problem is Solved
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="one_over_f_filter__01.png" alt=""> 
</div>
</div>
</body>
</html>